The potential mechanism of Guizhi Fuling Wan effect in the treatment of cervical squamous cell carcinoma: A bioinformatics analysis investigation

As a global malignancy with high mortality rate, targeted drug development for Uterine Cervical Neoplasms is an important direction. The traditional formula Guizhi Fuling Wan (GFW) is widely used in gynecological diseases. However, its potential mechanism of action remains to be discovered. We retrieved GFW and cervical squamous cell carcinoma (CSCC) targets from public databases. The protein–protein interaction network was obtained by string computational analysis and imported Cytoscape_v3.9.0 to obtain the core network and the top 10 Hub genes. Gene Ontology and Kyoto Encyclopedia of Genes and Genomes were used for enrichment analysis of the core network, and then molecular docking to verify whether the selected signaling pathway binds well to the core node. Finally, clinical prognostic analysis and expression differences of Hub genes were validated using the Cancer Genome Atlas database and R language. Our search yielded 152 common targets for GFW and CSCC. The interleukin-17 signaling pathway, tumor necrosis factor signaling pathway, and Toll-like signaling pathway were then selected for further molecular docking from the hub genes enrichment analysis results, which showed good binding. Among the Hub genes, JUN, VEGFA, IL1B, and EGF had a poor prognosis for CSCC. In conclusion, this study illustrates that GFW can have adjuvant therapeutic effects on CSCC through multiple targets and multiple pathways, providing a basis for further research.


Introduction
Uterine cervical neoplasms, mostly caused by persistent infection of the genital tract with 15 high-risk human papillomavirus types, are the most common malignant neoplasms of the female genital tract. [1,2]According to data provided by the International Agency for Research on Cancer in 2020, cervical cancer is the fourth most common cancer in women, after breast, colon and lung cancer, and is the fourth leading cause of cancerous death in women. [3]Uterine cervical neoplasms are found in the junction of old and new squamous epithelium and columnar epithelium of the uterine cervix, and squamous carcinoma is the most common tumor type, followed by adenocarcinoma and adenosquamous carcinoma. [4,5]urrently, the treatment plan for cervical squamous cell carcinoma (CSCC) is designed mainly based on the tumor stage, with neoadjuvant chemotherapy, radiotherapy and surgical resection complementing each other.Research on targeted drugs has also been the focus of CSCC treatment, and with the development of biological disciplines such as genomics, trials of the combination of multiple targeted drugs tend to transform from single target targets to multi-target networks. [6]n clinical practice, platinum drugs are important adjuvants, and the combination of carboplatin, paclitaxel, and bevacizumab has been shown to be effective. [7,8]This shows that the combination of drugs is a wise disease treatment option and research direction. [9]Chinese medicine has important research value in clinical practice as a coordinated overall multi-pathway and multi-target drug intervention.In recent years, traditional Chinese medicine has gradually been shown to exert significant antitumor and adjuvant effects in the treatment of cervical cancer by inhibiting cancer cell activity, regulating apoptosis, and modulating cell signaling variant pathways. [10,11]uizhi Fuling Wan (GFW), which is derived from the next volume of Jin Kui Yao by Zhang Zhong Jing, is a herbal formula composed of Cinnamomi Ramulus, Poria Cocos(Schw.)Wolf., Persicae Semen, Paeoniae Radix Alba, and Peony. [12][15] After animal model tests as well as clinical observation and analysis, GFW was found to improve the effects of uterine fibroids. [16,17]Meanwhile, it has been shown that GFW can reduce autophagy in granulosa cells of rats with polycystic ovary syndrome by restoring the PI3K/AKT/mTOR signaling pathway. [12]GFW was also shown to inhibit granulosa cell autophagy in polycystic ovary syndrome via H19/miR-29b-3p. [18]In addition to this, it has now been shown that extracts of GFW might be able to induce inhibition of target protein expression and restore sensitivity in cisplatin-resistant patients. [19]Therefore, GFW has shown its value in oncology treatment.Based on the existing studies of GFW for gynecological diseases and tumors, the present study tried to investigate the potential of GFW for the treatment of CSCC.However, as a compound, GFW has a complex system of chemically active components and potential targets of action, and it is difficult to elaborate its mechanism of disease treatment. [20]etwork pharmacology is a new field at the intersection of artificial intelligence and medicine, which is based on a large number of medical databases and computer algorithms to construct network relationship diagrams to illustrate the complex mechanisms of drug action. [21,22]The multi-target action in network pharmacology is consistent with the complex mechanism of action of diseases and herbal medicines, [23] which is free from the limitation of "one disease-one target-one drug" [24] and provides new ideas for the study of herbal formulations. [25]In order  to extensively verify the pharmacological mechanism of GFW action on CSCC, we will adopt the analytical method of network pharmacology. [26]The detailed workflow diagram of the study is presented in Figure 1.

Establish the databases of GFW target genes and CSCC-related genes
The relevant data on herbal medicines and diseases were obtained through the Traditional Chinese Medicine Systems Pharmacology (TCMSP) database (https://old.tcmsp-e.com/tcmsp.php)and the GeneCards database (https://www.genecards.org/).First, the TCMSP database was applied to search the relevant active ingredients of the 5 Chinese herbal medicines in GFW, and the screening criteria of oral bioavailability ≥30% and drug-likeness ≥0.18 were used to collect the potential action targets of each active ingredient respectively.Then, GeneCards was used to collect CSCC-related genes, and "Cervical Squamous Cell Carcinoma" was used as a keyword to search.

Constructing the relationship between active ingredients of Chinese medicine and diseases
The relevant drug targets collected from TCMSP were imported into the UniProt database (https://www.uniprot.org/)and converted to a common format to align them with the disease targets.GFW and CSCC common genes were obtained by Draw Venn Diagram (https://bioinformatics.psb.ugent.be/webtools/Venn/).

Construction of protein-protein interaction network
Import the GFW and CSCC common genes into the String database (https://cn.string-db.org/),select the "Homo sapiens" species, obtain the protein-protein interaction (PPI) network, and export it in " tsv " format file.Afterwards, the files were imported into Cytoscape _ v3.9.0 for visualization and analysis to generate PPI networks.

Acquisition of core networks
The network with the highest score was selected as the core network by analyzing the PPI network with the plug-in MCODE in  Cytoscape _ v3.9.0 software.The core network was then analyzed using CytoHubba, a plug-in in Cytoscape _ v3.9.0, to obtain the top 10 ranked Hub genes for follow-up analysis and research.

Gene Ontology and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis
The HIPLOT (https://HIPLOT.com.cn) is widely used for bioinformatics analysis as a powerful online data analysis website. [27]e used its Gene Ontology (GO)/Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis module to investigate the core network.The core genes were imported into the KEGG database in HIPLOT, and the species was selected as "Homo sapiens" with the P threshold set to 0.05.Then the R package "clusterProfiler" was used to enrich the genes to obtain valid information and obtain the results of GO and KEGG analysis.

Molecular docking
We selected from the relevant important signaling pathways enriched by KEGG that are closely related to the disease for the study.The targets of related signaling pathways were collected and collated, and gene targets with links to two or more signaling pathways were selected for follow-up studies.At the same time, the active ingredients represented by the targets were identified according to the previously established database and finally the Chinese herbal medicine to which the active ingredients belonged.The obtained results were imported into Cytoscape _ v3.9.0 software for visual analysis and molecular docking validation.Molecular docking requires pretreatment of target proteins and small molecules.The protein molecules obtained from the PDB database were modified by PyMOL software.Firstly, we remove the water of crystallization from the protein structure and use the "GetBox Plugin" plug-in to select the active pocket coordinates when the ligand is acquired for docking, then after removing the ligand, we obtain the original ligand of the target protein and export the protein molecule in the pdb.format.
Using AutoDock Tools to add polar hydrogen atoms and charges to proteins and to define atom types.The root of the ligand is determined by the "Torsion Tree" module and a twistable bond suitable for the ligand is selected.After completion, target proteins and small molecules in pdbqt.format will be exported for the next step of molecular docking.The obtained target proteins and small molecules in pdbqt.format will be imported into AutoDock Vina software for molecular docking.To obtain the best conformation of the complex, each molecule will be docked from 10 conformational angles.Finally, Discovery Studio 2020 software was used to perform the constitutive relationship analysis and visualize the results.

Prognostic analysis
We performed a relevant reanalysis and validation of the previously obtained top 10 Hub genes rankings.RNAseq data (level3) and corresponding clinical information of CSCC were obtained from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/).The log rank was used to test the Kaplan-Meier survival analysis, comparing the survival differences between the above two or more groups, and a timeROC analysis was performed to compare the predictive accuracy of the top 10 ranked Hub genes.For Kaplan-Meier curves, P values and hazard ratios with 95% confidence intervals were derived by logrank test and univariate Cox regression.All the above analysis methods and R packages were performed using R software version v4.0.3 (R Foundation for Statistical Computing, 2020).P < .05 was considered to be statistically significant.

Validation of differences in expression regarding Hub genes
RNAseq data (level3) and corresponding clinical information for 253 CSCCs were obtained from TCGA dataset (https://portal.gdc.com).We used R software v4.0.3 for statistical analysis and ggplot2 (v3.3.2) to visualize the expression differences of the 10 Hub genes obtained from the analysis.P < .05 was considered to be statistically significant.

Acquisition of GFW active ingredients and gene targets
The active ingredients of CR, PS, PRA, CM and PC were screened in TCMSP database with oral bioavailability value  ≥30% and drug-likeness ≥0.18,The results showed that there were 58 active ingredients and 185 gene targets in the 5 herbal medicines.The relevant active ingredients are shown in Table 1.The established GFW database was imported into Cytoscape _ v3.9.0 for visualization and analysis, and the active ingredient -gene target relationships were obtained as shown in Figure 2.

Common genes of GFW and CSCC
A total of 4533 targets related to CSCC were retrieved using GeneCards.The collected GFW database gene targets were changed into a common format using the UniProt database and then imported into Draw Venn Diagram together with CSCC genes to construct the Venn diagram, and 152 GFW and CSCC common genes were obtained, as shown in Figure 3.

Acquisition of PPI and core networks
As shown in Figure 4A, there are 150 targets and 2770 edges in the PPI network.The PPI network was then further analyzed using the MCODE plug-in, and the network with a score of 45.927 and the highest ranking was selected as the core network, as shown in Figure 4B.This core network has 56 nodes and 1263 edges.The core network was then processed with CytoHubba to obtain the top 10 ranked Hub genes with the highest association in the core network (JUN, VEGFA, IL6, CASP3, TNF, PTGS2, AKT1, IL1B, MMP9, and EGF), which is displayed in Figure 4C.

Results of enrichment analysis of GO and KEGG pathways
GO and KEGG enrichment analysis of 56 core targets was performed using the HIPLOT online analysis website to explore the mechanistic role of GFW in CSCC treatment.The results of GO analysis are shown in Figure 5A-C, involving biological process, cell composition, and molecular function, respectively.The results of the analysis of the top 15 pathways of KEGG are shown in Figure 5D, where we found that GFW treatment of CSCC was closely associated with IL-17 signaling pathway, tumor necrosis factor (TNF) signaling pathway and Toll-like receptor signaling pathways were closely related.

Molecular docking validation
Based on the P-value and the requirements of the study, we selected IL-17 signaling pathway, TNF signaling pathway and Toll-like receptor signaling pathway for the study.We imported a total of 28 targets from the 3 pathways into Cytoscape _ v3.8.2 for visualization and analysis, and the results are shown in Figure 6.
The binding energies of all the targets in the 3 pathways were ranked with the corresponding components, and the results are shown in Table 2.The top 5 combinations ranked from low to high were verified by molecular docking with the software AutoDock Tools, and the visualization of the best conformational fit obtained is presented in Figure 7.

The results of clinical prognostic analysis and expression difference validation regarding Hub genes
We used the TCGA database to reanalyze the Hub genes and the Kaplan-Meier survival curves obtained are shown in Figure 8. JUN, VEGFA, IL1B, EGF the 4 Hub genes had statistically significant differences in prognostic values, which was P < .05.For all 4 genes their hazard ratios was greater than 1, which means that there was a significant negative effect on prognosis.
The differences in gene expression between patients with CSCC and normal subjects are shown in Figure 9.We found significantly lower expression of JUN, VEGFA and AKT1 genes (P < .05)and elevated expression of GASP3, TNF, PTGS2, IL1B and MMP9 genes (P < .05),while there was no remarkable difference in IL6 and EGF gene expression changes.

Discussion
First, based on the ADME model, we screened 185 potential GFW targets through the database, and obtained the PPI network with 152 common genes after taking the intersection with the valid targets of CSCC.Then, after obtaining the core network and the top 10 ranked Hub gene, we performed genetic differential analysis and prognostic analysis to verify the correlation between the obtained Hub gene and the disease.
Interestingly, based on the available clinical data, it can be found that the expression of JUN, VEGFA, AKT1, GASP3, TNF, PTGS2, IL1B and MMP9 are significantly different from normal subjects, while the prognosis of JUN, VEGFA, IL1B and EGF is significantly different.For JUN, a proto-oncogene, it is highly susceptible to loss and translocation malignant tumor development.The expression product of JUN is involved in the composition of transcription factor AP-1, which is associated with important functions such as epithelial cell proliferation, apoptosis and histomorphogenesis. [28,29]The expression of JUN is significantly downregulated in CSCC compared to normal subjects, and the prognosis of the high JUN expression group in CSCC prognostic analysis is poorer, and the expression products of JUN have an influence on the developmental process of CSCC.VEGFA induces tumor angiogenesis by regulating the expression of a series of downstream products to achieve enhanced vascular endothelial cell viability and the formation of an abnormal tumor vascular system.[32] IL1B, a member of the interleukin family, is a key pro-inflammatory factor involved in cell proliferation, differentiation and apoptosis, which can promote tumor growth and metastasis.High expression of IL1B can lead to poor prognosis in CSCC. [33,34]As for EGFR, it has been an important target for cancer therapy, [35] for which EGF is a ligand.EGF/EGFR has an important role in promoting cell division and enhancing epidermal cell viability.In vitro experiments demonstrated that the growth of Hela cervical cancer cells remained dependent on high EGFR expression. [36]It has been shown that upregulation of EGF expression plays a role in the development of CSCC. [37]he expression of EGF is higher in CSCC than in the general population, and the prognosis of the high expression group is significantly worse.Due to the lack of clinical data, the relevance of some Hub gene to CSCC remains to be explored.The above findings suggest that the potential core targets of GFW action are closely associated with CSCC, and the efficacy of GFW and CSCC has some research value.The individual targets are interconnected in the core protein network and often act by regulating multiple signaling pathways.Next, we performed KEGG enrichment using core genes.From the top 15 pathways with the highest enrichment scores, IL-17 signaling pathway, TNF signaling pathway, and Toll-like signaling pathway were selected for this study.IL-17 family includes 6 ligand members and 5 receptor members from IL-17A to IL-17F, and the specific proximal adaptor downstream of IL-17R is Act1, which activates downstream pathways including NF-κB, TNF, and Toll through a series of phosphorylation and ubiquitination, activating downstream pathways including NF-κB, MAPKs and C/EBPs [38][39][40] to regulate cellular activity.Notably, as a classical pro-inflammatory pathway, the IL-17 signaling pathway is also closely associated with tumors, where it can mediate immunosuppression and angiogenesis.A previous study showed that HPV16 oncoprotein can promote angiogenesis and CSCC development by inducing the production of IL-17A+γδ T cells. [41]The TNF signaling pathway cascades with multiple pathways and plays an important role in pathophysiological processes such as cell proliferation, apoptosis, and pro-inflammation. [42,43]Currently, it has been that TNF-α can inhibit normal cervical epithelial cells, and stimulate the proliferation of malignant cervical epithelial cells, providing a growth advantage for CSCC. [44]The NF-κB signaling pathway acts as a downstream activation pathway of the TNF signaling pathway, [45,46] where NF-κB has been shown to be constitutively activated during the development of CSCC. [47]The TLR is involved in important molecules in innate immune defense and are thought to have an important role in the clearance of HPV virus infection. [48]In CSCC-related studies, TLR4 upregulation and TLR2, TLR7 downregulation were found to be associated with poor prognosis in CSCC, [49] providing an idea for TLR agonist studies.This can prove that the key signaling pathways enriched in core targets are closely associated with CSCC, and the mechanism of action of GFW has been initially revealed.
Finally, to further validate, we selected 28 core targets of IL-17 signaling pathway, TNF signaling pathway and Toll-like signaling pathway, and molecularly docked them with the corresponding active ingredients of GFW.From the results (Table 2), quercetin, kaempferol, taxifolin and other active ingredients showed high binding efficacy to the relevant targets in the signaling pathway, providing strong evidence for the elucidation of the mechanism of action of GFW on CSCC in terms of molecular docking.Notably, quercetin is an important molecule with anti-tumor effects, [50] and in vitro experiments have demonstrated the potential of quercetin to enhance the sensitivity of human cervical cancer cells to cisplatin-induced apoptosis, [51] playing a potential adjuvant therapeutic role.By establishing a "signaling pathway-active component-Chinese medicine" network, we can find that the 5 major components of GFW, PRA, CR, PS, CM, and PC, all play important roles.
In summary, this study provides a preliminary description of the mechanism of action of GFW on CSCC by establishing a closely linked "active ingredient-target-signaling pathway" network through network pharmacology and bioinformatics analysis.The study demonstrated the potential value of GFW as an adjuvant treatment for CSCC and provided an idea for further research in the future.Of course, this study is based on the existing database and previous studies, and the conclusions obtained lack experimental validation and need to be further explored in the future.

Conclusion
In this study, we used network pharmacology and molecular docking techniques to verify the specific mechanism of GFW action on CSCC with the concept of "multi-target-multipathway."The results demonstrated that GFW can act on the relevant targets through IL-17 signaling pathway, TNF signaling pathway and Toll-like signaling pathway, and the relevant active ingredients bind well to the important targets in the pathway.This result reveals that GFW has potential adjuvant therapeutic effects on CSCC and provides a theoretical basis for further experimental validation.
F), AKT1 (G), IL1B (H), MMP9 (I), and EGF (J), in that order.The HR (High exp) represents the risk coefficient of the high expression group relative to the samples in the low expression group, if HR > 1 means the gene is a risk factor (the higher the expression, the worse the prognosis), if HR < 1 means the gene is a protective factor (the higher the expression, the better the prognosis); 95% CL represents the HR confidence interval; Median time represents the time (i.e., median survival time) corresponding to the survival rate at 50% in both the high expression and low expression groups in years.HR = hazard ratios.

Figure 1 .
Figure 1.Network pharmacology and bioinformatics workflow diagram.

Figure 3 .
Figure 3. Venn diagram about the common goals of GFW and CSCC.There are 152 common goals between them.CSCC = cervical squamous cell carcinoma, GFW = Guizhi Fuling Wan.

Figure 2 .
Figure 2. Network diagram of active ingredients of GFW in relation to gene targets.Green circular nodes represent gene targets; blue diamond-shaped nodes represent corresponding active ingredients.GFW = Guizhi Fuling Wan.

Figure 4 .
Figure 4.The analysis of the common goal relationship between GFW and CSCC.(A) PPI network diagram.(B) Core network: it consists of a combination of the top 50 common targets, in which the darker the node color, the higher the ranking.(C) The top 10 ranked Hub genes in the core network.Nine red circular nodes represent the highest ranked genes, they are all the top-ranked genes, and the yellow circular nodes in the middle represent the lowest ranked genes.CSCC = cervical squamous cell carcinoma, GFW = Guizhi Fuling Wan, PPI = protein-protein interaction.

Figure 5 .
Figure 5. Enrichment analysis of GO and KEGG pathways.(A-C) Results of BP, MF, and CC analyses obtained by GO enrichment, respectively.(D) are the results obtained from KEGG enrichment analysis.In the bar or bubble plots, the smaller the P-value, the higher the enrichment c-degree, while the longer the bar or the larger the bubble, the more the number of enriched genes.BP = biological processes, CC = cell composition, GO = gene ontology, KEGG = Kyoto Encyclopedia of Genes and Genomes, MF = molecular functions.

Figure 6 .
Figure 6.Signaling pathway, target, active ingredients and herbal medicine relationship diagram.(A) Diagram of the relationship between signaling pathways and corresponding targets.Green triangular nodes represent signaling pathways; purple circular nodes represent gene targets that are linked to one of the signaling pathways; blue circular nodes represent gene targets that are linked to two of the signaling pathways; yellow circular nodes represent gene targets that are linked to all 3 signaling pathways.(B) Relationship diagram between gene targets and corresponding active ingredients, orange rectangular nodes represent active ingredients.(C) Relationship diagram between the active ingredient and herbal medicine, dark blue rectangular nodes represent herbal medicine.

Figure 8 .
Figure 8. Plot of prognostic analysis of the first 10 Hub genes.The relevant analyses were performed for JUN (A), VEGFA (B), IL6 (C), CASP3 (D), TNF (E), PTGS2 (F), AKT1 (G), IL1B (H), MMP9 (I), and EGF (J), in that order.The HR (High exp) represents the risk coefficient of the high expression group relative to the samples in the low expression group, if HR > 1 means the gene is a risk factor (the higher the expression, the worse the prognosis), if HR < 1 means the gene is a protective factor (the higher the expression, the better the prognosis); 95% CL represents the HR confidence interval; Median time represents the time (i.e., median survival time) corresponding to the survival rate at 50% in both the high expression and low expression groups in years.HR = hazard ratios.

Table 1
Potentially active compounds in Guizhi Fuling Wan.

Table 2
The combination of the best docking model energy.